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We develop the relativistic theory of hydrodynamic fluctuations for application 
to high energy heavy ion collisions. In particular, we investigate their effect on the 
expanding boost-invariant (Bjorken) solution of the hydrodynamic equations. We 
discover that correlations over a long rapidity range are induced by the propagation 
of the sound modes. Due to the expansion, the dispersion law for these modes is 
non- linear and attenuated even in the limit of zero viscosity. As a result, there is a 
non-dissipative wake behind the sound front which is generated by any instantaneous 
point-like fluctuation. We evaluate the two-particle correlators using the initial con- 
ditions and hydrodynamic parameters relevant for heavy-ion collisions at RHIC and 
LHC. In principle these correlators can be used to obtain information about the 
viscosities because the magnitudes of the fluctuations are directly proportional to 
them. 

I. INTRODUCTION 

The success of relativistic hydrodynamics in describing the fireball created in ultrarela- 
tivistic heavy ion collisions opened the possibility to study the properties of strongly inter- 
acting matter at extremely high temperatures and densities near thermal equilibrium. We 
know from lattice simulations of quantum chromodynamics (QCD) that strongly interacting 
matter at temperatures above the crossover at Tc ~ 165 MeV is a quark-gluon plasma [III2]- 
Lattice QCD is able to predict stationary thermodynamic properties of the quark-gluon 
plasma, such as the equation of state, but is presently unable to make reliable predictions 
for dynamical properties, such as transport coefficients. 
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A remarkably small value of the shear viscosity t] in the natural units of the entropy 
density s, Atttj/s < 2.5, has been deduced from comparison of the results of relativistic 
viscous fluid dynamics simulations with data from Au+Au collisions at the Relativistic 
Heavy Ion Collider (RHIC) [3J. This result is interesting because it is smaller by at least a 
factor 5 than the value oirj/s calculated in thermal perturbation theory at leading order |1] 
and not far away from the value found earlier in a large class of strongly coupled non-abelian 
gauge theories [5] . It would thus be desirable to confirm the inferred experimental result for 
the shear viscosity by other methods. 

Due to the fluctuation-dissipation theorem, the shear and bulk viscosities not only control 
the dissipative properties of a fluid in the limit of small velocity gradients, but they also 
control the magnitude of hydrodynamic fluctuations in the fluid. ^ Thus it is interesting 
to explore whether fluctuations in the density and flow velocity of the fluid can be used 
to deduce the value of the shear viscosity from experimental data. The purpose of our 
work is to lay the foundations for quantitative investigations of this idea by identifying the 
sensitivity of correlation observables to the hydrodynamic fluctuations. 

Before proceeding, let us identify four major sources of density fluctuations in relativistic 
heavy ion collisions: 

(a) Initial state fluctuations: These are the result of quantum fluctuations in the densi- 
ties of the two colliding nuclei and fluctuations of the energy deposition mechanism. 
They appear as event-by-event fluctuations in the energy density and flow velocity 
distributions at the onset of the hydrodynamic regime. These fluctuations and their 
phenomenological ramifications have recently been studied extensively [7H20] because 
they may be responsible [2T] for the angular correlations of particle emission observed 
in the heavy-ion experiments [221427] . The power spectrum of the final-state angu- 
lar correlations induced by initial-state fluctuations may provide information about 
the speed of sound and the shear viscosity of the matter produced in the heavy-ion 
collision 128] . Longitudinal fluctuations and correlations among the initial-state 
angular fluctuations have been investigated [T71 122] ■ The initial state correlations 
over large rapidity intervals have been subject to studies in connection with the ridge 

^ This fact is also represented by Kubo formulas, relating viscosities to correlators of stress-energy tensor, 
and underlies the approach taken in Ref. [5] to study bulk viscosity. 
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phenomenon j30l439] . 

(b) Hydrodynamic fluctuations: These are the result of finite particle number effects in 
a given fiuid cell. This leads to local thermal fiuctuations of the energy density and 
fiow velocity which propagate according to the hydrodynamical equations. According 
to the general theory of hydrodynamical fiuctuations [40], the squared amplitude of 
these fiuctuations is proportional to the viscosity. These fiuctuations are the focus of 
our paper. 

(c) Fluctuations induced by hard processes: Energetic partons, which have been scattered 
in the initial collision of the two nuclei, can propagate through the quark-gluon plasma 
where they lose energy. To the extent that this energy is thermalized, it acts as a source 
term for the hydrodynamical equations. The space-time shape of this source term has 
been calculated in the weak and strong coupling limit . If the shear viscosity of 
the plasma is as low as inferred from the RHIC data, these sources will excite Mach- 
cone shaped perturbations in the expanding fiuid [SI |15]. It is presently not clear 
whether these lead to observable phenomena after freeze-out ^6] . 

(d) Freeze-out fluctuations: Event-by-event fiuctuations may also be caused by finite par- 
ticle number effects during and after the freeze-out of the hydrodynamically expanding 
fluid. 

The main purpose of this article is to develop and apply the relativistic theory of hy- 
drodynamical fluctuations to the evolution of the quark-gluon plasma formed in relativistic 
heavy-ion collisions. Although the relativistic generalization of the text-book |1D] theory of 
the hydrodynamic fluctuations have been considered before in a different context [17], for 
completeness we outline the derivation in Sec. [Tlj We then apply the resulting stochastic 
hydrodynamic equations to the simplest example of a boost-invariant (Bjorken) flow. Our 
purpose is to illustrate the application of the theory in a most transparent, yet phenomeno- 
logically meaningful setting. 

We are able to obtain a number of closed form analytic results which demonstrate im- 
portant phenomenological consequences of the theory. In particular, we flnd that hydro- 
dynamic fluctuations during the early phase of the expansion naturally induce correlations 
across large rapidity intervals. It is usually assumed that such correlations, observed in 
experiments, must be produced before equilibration. 
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A full implementation of the framework we develop here will eventually require numerical 
solution of stochastic hydrodynamic equations in three spatial dimensions, which we defer 
to future studies. 

II. THEORY OF HYDRODYNAMIC FLUCTUATIONS 

Hydrodynamics is an effective theory that describes the long wavelength and low fre- 
quency space-time evolution of the densities of a few conserved quantities such as energy, 
momentum, electric charge and baryon number. In the case of spontaneous breaking of con- 
tinuous symmetries it also describes the evolution of the phases of order parameters. These 
hydrodynamic variables are defined as average values of the corresponding local, space-time 
dependent, coarse-grained operators. The coarse-grained averaging is performed over dis- 
tances and times that are small compared to the macroscopic scales of interest but large 
compared to the microscopic scales, such as mean free paths and mean free times between 
collisions. The hydrodynamic variables evolve according to the deterministic equations which 
follow from the conservation equations obeyed by the corresponding operators. 

Fluctuations and correlations in the hydrodynamic variables can be characterized by 
averaged values of the products of the operators at different space-time points. Although the 
fluctuations themselves occur on microscopically short space-time scales, these fluctuations 
are correlated not only on short space-time scales but also on macroscopically large scales. 
This can be understood as a result of the diffusion or propagation of each fluctuation at any 
earlier time to later times. Such propagation over long times and distances, and thus the 
long-range behavior of correlation functions, is described by hydrodynamics. 

A. Hydrodynamic variables and equations 

We consider a general case of a system with 5 conserved quantities: energy, charge and 
three momentum components. In the case of QCD we can think of the charge being the 
baryon number. The five equations are the conservation equations for the energy-momentum 
tensor, d^Ti^^ — 0, and the current conservation equation d^J'^ — 0. 

The energy-momentum tensor and current densities for a fluid in thermal, chemical and 
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mechanical equilibrium are 

7]^^^^ = -Pg'"' + wu'^u''; = hbu^" (equilibrium) (1) 

Here P is the equilibrium pressure at given energy density e and baryon density ns, w = 
P + e = Ts + fisnB is the local enthalpy density, is the baryon chemical potential and 
is the local flow 4-velocity. The metric Qf^u is (+, — , — , — )• 

The non-equilibrium corrections to these expressions, AT^'^ and AJg, are proportional, 
at lowest order, to flrst derivatives of the local quantities with coeflicients given by the shear 
viscosity t], bulk viscosity and thermal conductivity x- Explicit expressions may be found 
in textbooks [IHl Il9] . Local thermal fluctuations are described by the additional stochastic 
terms S^"" and I^"". 

= riBu" + AJ^ + (2) 

In the following we shall determine the correlation functions of the stochastic terms. Since 
the source of the fluctuations is local, these correlation functions are delta-functions in 
space and time. The amphtude of these source terms is flxed by the fluctuation-dissipation 
theorem. 

In practice, the idea is to consider the stochastic terms as given functions of space and 
time and to solve the fluid equations of motion to flrst order in it. Quantities which are linear 
in the S'^" will average to zero, where the average is taken over the ensemble of fluctuations. 
Quantities which are quadratic in the S'^'^ may have non-zero average values which we must 
determine. 

The form of the hydrodynamic equations depend on the definition of the local flow velocity 
u'^. There are two common choices; we discuss each in turn, including the modifications 
necessary to incorporate fiuctuations. 

B. Eckart approach 

The Eckart approach is a convenient choice if we want to compare with the non-relativistic 
limit. In this approach is the velocity of baryon number fiow. The dissipative terms must 
satisfy the conditions A Jg = and u^UyAT'^^ = 0, the latter following from the requirement 



that 7"°° be the energy density in the local (baryon) rest frame. The most general form of 
AT^"' is 

AT^- = ATZ + ATZt (3) 

where 

ATZ = V i^'W" + A'^u^) + (§77 - C) h^"" (d ■ u) (4) 
is the viscous part and 

ATTat - X {h^'-u^ + hT'-u^) (5) 
is the heat conduction part. Here 

hi''' = ui'u'' - gi"" (6) 

is a projection tensor normal to u^^, 

A^ = d^-u^{u-d) (7) 

is a derivative normal to m^, and 

= -d^T + T{u-d)ua (8) 

is a four-vector whose nonrelativistic limit is q = VT". The entropy current is 

= sui" + ^u^ATi"" . (9) 

By using energy-momentum conservation d^T^^'^ — the divergence of the entropy current 
can be put in the compact form 

d^s'^ = AT'^^a^ . (10) 
For some purposes it is better to express this divergence as 



d^s^ = ATI"' 



(11) 



Substituting the explicit form of AT^'' into the above expression gives 

d^si' = ^ [{Ai^u" + A^ui") + {d ■ u)] ' 



+ ^{d-uf + ^^hi^^q^q. . (12) 
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In the local rest frame this is 

+ ^{V-uf + ^^{VT + Tuf . (13) 

The term Tii is a relativistic correction to VT, being smaller by a factor of in physical 
units. All three dissipation coefficients must be non-negative to insure that entropy can 
never decrease. 

It is useful to decompose S^^ into a piece associated with viscosity and another piece 
associated with heat conduction. Overall we must require that u^UuS^" = 0, just like AT'^^, 
so that in the local rest frame of the fluid equals e = Tj^^j^j. Then if we are given S^" 
with this property we can define 

-^he^at = S'-'^^'aU" + S^^^M^M^ (14) 

and 

This decomposition is unique. 

We follow Section 88 of [IQ] on hydrodynamic fiuctuations. In the general theory of 
quasi-stationary fiuctuations, presented in one considers the set of equations 



^a = - ^ labXb + Va (16) 
b 



which gives the response of the set of variables Xa to the driving terms Xa and to the ya, 
which represent random fiuctuations. The time rate of change of the entropy is 

S=-Y,XaXa. (17) 

a 

In order for the probability distribution of fiuctuating variables to agree with the thermo- 
dynamic distribution given by e"^, the noise autocorrelations must be given by 

{VaitMh)) = (7afe + lba)5{h ' ^2) ■ (18) 

This general framework needs to be applied to the present situation. 
The time rate of change of the total entropy of the system is 



dS_ 
~dt 



(19) 
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Coarse graining is performed in the usual way with cell volumes AV. Since viscosity and 
heat conduction are independent physical processes, it is natural to make the identifications 

i2 ^ ATZ, . (20) 
Comparing to the rate of entropy change allows us to deduce that 
Xi -^[A^u, + A,u^]AV, 



AV. (21) 



Next the Onsager coefficients '^ab can be determined. 

711 = 2T [rjh^-h''^ + i (C - i^) h^'h"^] ^ , 

722 = 2xT^ [h^^u^u^ + K'^u^u''] ^ . (22) 

The 7ii is made unique by the requirement that it vanish when any of its indices is contracted 
with the four-velocity. The 712 and 721 are zero as expected. 

Different coarse grained cells are independent. Then the factor 1/ AV goes over to a Dirac 
delta function in position space. The correlation functions are easily written down (after 
acknowledgement that they must have certain symmetries in the Lorentz indices) . They are 

+ (C-i7?)«"T(5(xi-X2) (23) 



and 



+ h^'^u'^u'' + K^'u^'u^] 5{xi - X2) (24) 



and 



('5v^:(^i)'5£t(^2)) = 0. (25) 

When the viscous correlation function is evaluated in the local rest frame it will vanish 
unless all of the indices are spatial. With ixu = ik and a(3 = Im we get 

{Sill^{x,)S!Zix2)) = 2T [rj (SuSkm + SimSki) + (C - h) S^kSlm\ S{x, - x^) (26) 
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which is exactly the expression in pU]. When the heat correlation function is evaluated in 
the local rest frame it will vanish unless each Sheat has one spatial and one temporal index. 
With /iz/ = Oi and a(3 = Oj we get 

(^h:at(^i)^hL(^2)) = 2xT%,S{x, - X,) (27) 

which also agrees with the corresponding expression in [30] • Since these correlation functions 
reduce to the known ones in the local rest frame, and since they are constructed from tensors, 
they are obviously valid in any frame of reference. 



C. Landau-Lifshitz approach 

The Landau-Lifshitz approach is the most convenient and frequently used approach for 
ultrarelativistic heavy ion collisions. In this approach is the velocity of energy transport. 
The dissipative part of the energy-momentum tensor satisfies m^AT'^'" = 0, and AJg is not 
constrained to be zero. In this case the most general form of the energy-momentum tensor 
is 

AT^- = AT4" = 7] {A^u' + A'^M^) + - C) h^'^d ■ u . (28) 
The baryon current is modified by 

A = orTA^ (P^xb) , (29) 

where a is the (baryon) charge conductivity. The modification to the current satisfies 
Ufj_AJ^ = 0. This means that n-g is the baryon density in the local rest frame. 
The entropy current in this approach is different, being 

s^' = su" - /3/iBA . (30) 

Using baryon number conservation, 9^Jg = 0, we can write 

d^s" = {su") + /3/iB9^ {riBu") - AJ^d^ {(3^ij,) . (31) 

By using energy-momentum conservation in the form u^d^T"^ = 0, this can be written in a 
way convenient for future use. 

" 1 



d,s" = ATZ 



2T 



{A^Uy + A^u^) 



+ A [VA'^ (/3/iB)] . (32) 



10 



Compared to the Eckart frame there is no change in the viscous part associated with shear 
and bulk viscosities. Therefore it can again be written in the symmetric form 

+ aW'^A^(^//B)A,(^//B) . (33) 

The part due to charge conductivity seems to be different than the part due to heat 
conduction in the Eckart frame, but it is not. Using energy-momentum conservation in 
the form ^a^^i^^deai ~ 0, which is vahd to zeroth order in the dissipative coefficients and 
sufficient for this purpose, and dP — sdT + nsd/iB, one finds 

w 

Aa(/3/^B) = ^ga. (34) 

This can be inserted into the expression for the divergence of the entropy current to obtain 
exactly the same expression as in the Eckart frame, provided that the charge conductivity 

a is related to the heat conductivity x, by 

a = xT{nB/w)^ (35) 

which corresponds to the Franz- Wiedemann law. 

The fluctuations S'^ = and must satisfy the conditions u^j^S'^^ = and Ufj^I^ — 
for the reasons mentioned above. 

The time rate of change of the total entropy of the system is 

^ = Jd'x J^ATZ ^ {A,u. + A^u,) + A [h,.A-^ (/3/.b)] } • (36) 
It is natural to make the identifications 

X2 AJ^. (37) 
Comparing to the rate of entropy change allows us to deduce that 

Xi -^[A^u, + A,u^]AV , 

X2 ^ - VA" (/9/^b) Al^ . (38) 
Next the 7a(, can be determined. 

711 = 2T [Tjh'^'^h''^ + i (C - Iv) h^'h'^^] ^ , 

722 = <7W^^. (39) 
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The 7ii is made unique by the requirement that it vanish when any of its indices is contracted 
with the four-velocity. The 722 is made unique by the requirement that w^AJg = 0. The 
7i2 and 721 are again zero as expected. 

The correlation function for the viscous part (5'^s(a;i)S'^f (X2)) is exactly the same as in 
the Eckart approach. The mixture (S'^g(a;i)J"(a;2)) is zero. The correlation function for the 
baryon current is 



When the baryon current correlation function is evaluated in the local rest frame it will 
vanish unless both indices are spatial. Then 



This completes the generalization of the theory of hydrodynamic fluctuations to the rela- 
tivistic domain. 

III. FLUCTUATIONS IN BOOST INVARIANT HYDRODYNAMICS 

In this section we consider, as an example, application of the stochastic hydrodynamic 
equations derived in the previous section to the hydrodynamic fluctuations around Bjorken's 
boost-invariant solution of relativistic hydrodynamics [51]. Unlike the thermal fluctuations 
around a stationary equilibrium solution, which are well-known, the correlations induced by 
hydrodynamic fluctuations on a non-stationary solution have not been discussed in the liter- 
ature to our knowledge. Although this example is not entirely realistic or directly applicable 
to data, it is semi-analytic in nature and allows us to gain experience with these fluctuations 
and with the issues that may arise in more realistic, multi-dimensional calculations. 

Here we shall consider only fluctuations of temperature (or energy density) and flow ve- 
locity and neglect the effects of the baryon number fluctuations. For highly relativistic heavy 
ion collisions at LHC and the top range of RHIC energies this is a reasonable approxima- 
tion because the smallness of the baryon chemical potential /iB suppresses mixing between 
baryon charge and energy density fluctuations. 

In this example we shall focus on longitudinal flow fluctuations by integrating all densities 
over the coordinates x and y perpendicular to the beam or z axis. This effectively reduces 
the dimensionality of the problem to (1 + 1). Thus, our example is different from the 




(40) 




(41) 
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treatment in the existing literature in at least in two aspects: (i) we consider liydrodynamic 
fluctuations, not initial state fluctuations; and (ii) we consider longitudinal correlations, not 
azimuthal ones. We shall briefly discuss transverse correlations in Appendix |B] but defer 
their detailed study to further work. 

It is convenient to view the Bjorken boost-invariant flow in Bjorken coordinates: proper 
time T and spatial rapidity ^. 

r = y/t'^ — z'^ 
^ = tanh~\z/t) 
t = rcosh^ 

z = rsinh^ (42) 

The average values of hydrodynamic quantities depend only on r while fluctuations, after 
integration over the transverse coordinates x and y, depend on both, r and ^. The flow 
velocity is given by = x^/t + 6u^, where the last term denotes the fluctuations. We 
express the fluctuations of the longitudinal flow in terms of the rapidity variable u which 
we deflne as 

M° = cosh(^ + w(^,r)) 

= sinh(e + a;(e,r)). (43) 

The local pressure depends on the temperature which in turn depends on both coordinates. 
The average value of T depends only on the proper time, but fluctuations of T depend on 
both coordinates. Therefore 

T = To(r) + 5T(e,r) 
P = Po{t) + 6P{^,t) 

e = eo(r) + 5e(e,r), (44) 

where the subscript refers to the average value of the function. Obviously all variations 
are related to variations in the temperature on account of the equations of state. 

6e = cv{T)6T 

5s = '-^5T 
T 

5P = s{T)6T 

Sw = 6e + 6P. (45) 
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Here cv{T) = de/dT is the heat capacity per unit volume. For the case of zero chemical 
potentials, as we are considering here, cyiT) = s{T)/vs^(T). 

The noise term satisfies u^S^^ = and is symmetric in its indices. Due to the reduced 
(1+1) dimensionality of this model this condition allows us to express the noise in terms of 
a single scalar function / as 

S^^ = w{T)fi^,T)h^'' (46) 

where h^'^ was defined in Eq. (|6|; the factor of w{t) is included to make / dimensionless 
and to simplify subsequent formulas. For the same reason the viscous term can be expressed 
in terms of a single function 

ATZ = -{lv + C){d-u)h'^-^. (47) 

The fiuctuations mentioned above will be linear functionals of /. Their average values will 
be zero since (/) = 0. The fluctuations in those observables will be determined by 

(/(ei, ri)/(6, r,)) = -^^^^ [Wi] + C(ri)] S {n - r,) S (^i - 6) (48) 



on account of Eq. (23). In this case the delta-function in the transverse coordinates 
6 (x^i — x_L2) is replaced with 1/A, where A is the effective transverse area of the colliding 
nuclei (for noncentral collisions it would depend on the impact parameter in the usual way) . 

A. Hydrodynamic equations 

The hydrodynamic equations of motion can now be written out using any one of several 
standard methods. There are two independent scalar equations of motion each of which is 
first order in derivatives. In the absence of fiuctuations, one of them is satisfied automatically 
due to the assumption of boost invariance. When dissipation is neglected the nontrivial 
equation simply expresses entropy conservation for ideal fiuid fiow 

d(Ts) , 
^-^ = (49) 
dr 

and has the solution s(r) = s{to)to/t where tq is the initial proper time (when thermalization 
first is achieved). Once dissipation is included the equation becomes more complicated. 
Defining 

z/ = (47^/3 + OA, (50) 
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the equation is 



d{Ts) 

dr 



us 



(51) 



meaning that the entropy per unit rapidity interval, tsA, increases due to dissipation. The 
exphcit solution requires knowing the relationship between s and T, in other words the 
equation of state, plus the temperature dependence of u. For example, using an equation of 
state with v^^ = dP/de = constant, and with u = constant, the solution is 



T(t) 



v,^-llu 



(52) 



where 7^ = 1 / ^/\ 



Compared to the inviscid case the temperature decreases more 



slowly, assuming it starts from the same value. 

Now we account for fluctuations by adding noise. At this point, we make no assumption 
about the form of the equation of state or the temperature dependence of u. The two 
independent equations that follow are 



d5e _ „ 5{us) 

T— + Sw + Wf ^ + 



and 



r 



9_ 



u w 



us 

T 



dr 



+ 2uj (w 



T 

d_ 



us 



w 



T 



5P + wf 



5{us) 







us d'^uj 



(53) 



(54) 



In deriving these equations it is helpful to make use of Eqs. (46) and (47). On account of 



the reduced dimensionality, as reflected in (46) and (47), Eqs. (53) and (54) follow from the 
equations of motion of a perfect fluid (/ = and = 0) by the replacements 



P ^ P + wf 



us 

T 



1 + 



- 



while 6e is unchanged. The fluctuations 6P, 6s, 6e, and 6w can all be expressed in terms of 
a new dimensionless variable 

p = 6s/ s, (55) 

so that 6e = wp and 6P = Vg'^wp. Hence the pair of equations (53) and (54) will determine 
the two independent dimensionless variables p and u. 

Since the unperturbed solution is boost-invariant and independent of ^, it is advantageous 
to use the Fourier transform 



X{k,T) 



d^e-^'^Xi^,' 



(56) 
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for any variable X. Note that the wavenumber k is dimensionless. With this transformation 



Eqs. (53) and (54) become a pair of coupled first order linear differential equations. The 



solutions for the dimensionless variables can be expressed as 

p{k,T) = 

and 



dr' ~ 

— G,{k;T,T')f{k,T') 



TO 



ijj{k, r) 



dr' ~ 

-GUk;T,T')f{k,T') 



r 



(57) 



(58) 



Note that p(/c, tq) = and uj{k, tq) = so that there are no fluctuations in the initial 



conditions, although they could easily be incorporated (see Sec. IV C). 

The problem reduces to finding the (dimensionless) Green functions Gp and followed 
by quadrature. Averaging is performed by use of 

2z/(ri) 



(/(A;i,n)/(fc2,r2)) 



In fc-space the correlators are 



Atiw{ti) 



6{Ti-T2)2^6{ki + k2) . 



(59) 



X(A;i,ri)r(A;2,r2)) = ^5(A;i + fc2) 



min(Ti,r2) 



dT2u{r) 
w{r) 



Gxiki;n,T)GY{k2;T2,T) (60) 



TO 



where X and Y can be either p or u. The correlator in ,^-space is obtained by a Fourier 



transform (56) 



For the most part we shall be interested in the equal-(proper)time correlation function 
at the freeze-out time Tf, which can be written as 

Tf 

Cxy(ei-6;rf) = (X(ei,rf)F(6,rf)) = 4 / G'xy(ei - 6; n, r) , (61) 

A J T-^ W[T) 

TO 

where the Fourier transform of Gxy{^', ^f, t) is given by 



Gxrik; Tf, r) = Gx{k; Tf, T)GY{-k; Tf, r) . 



Thus 



Gxy{^i -6;n,r) 



rfeG'x(6-e;rf,r)Gy(6-e;rf,r) 



(62) 



(63) 



This equation shows directly that a fluctuation at point C, at time r induces a correlation 
between points and ^2 at later time T{ via a hydrodynamically propagating response given 



by Eqs. (57) and (58). 
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B. Inviscid case 

For the sake of clarity we shall first present the case where we neglect the contribu- 
tion of dissipation (viscosity) in the equations of motion and then later consider viscous 
corrections. Of course we cannot simply set the viscosity to zero since, according to the 



fluctuation-dissipation theorem expressed by Eq. (23), we would also eliminate fluctuations. 
However, if the viscosity is small and the flow is close to ideal, as is the case for heavy ion 
collisions, the contribution of viscous terms to the correlators is limited to the vicinity of 
singularities. These singularities correspond to unattenuated propagation of sound shocks 



which the viscosity will smear out, as we shall quantify in Sec. IIIF 



After Fourier transformation, Eqs. (53) and (54) become 



T^ + iku + f = (64) 
or 

+ (1 - ^s^) w + ikv^'^p + ikf = (65) 

There are at least two different methods to solve these equations, each having their own 
merits. They must, of course, yield the same solution. We outline each in turn. 

The pair of equations can be combined into a Langevin equation for the two-component 
vector 

^ ■ ^^^^ 

The Langevin equation takes the form 

+ D'4^ + h = 0, (67) 

or 

where the drift and noise terms are given by 

D = Do=i^^ I , h={' \f. (68) 




-n = 0, 

















The solution to these equations for arbitrary noise and drift, with the given initial conditions, 
can be written as 

^(fc,r) = - / —Uik;T,T')hik,T') (69) 
where U is the evolution operator satisfying 

^ dU{k^T,T') ^ ^^^^ ^ ^ ^^^^ 
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subject to the condition U{k;T,T) = 1. Explicitly 



-r J II 



(71) 



where T is the time ordering operator. Comparing Eq. (69) to Eqs. (57) and (58) allows 
for determination of Gp and G^, namely 



Guj{k; r, t') 



U{k;T,T') 



1 

ik 



(72) 



The second method is to eliminate one of the variables in favor of the other to arrive at 
a single second order differential equation. Elimination of Cj results in 



9r2 



dr 



0. 



(73) 



Denote the two independent solutions to the homogeneous equation, when / = 0, by pi 
and p2- The function Gp is constructed from a linear combination of the two homogeneous 
solutions as 

Gp{T, t') = ai(r')pi(r) + 52(r')p2(r) . (74) 



The functions di and 0.2 are determined by substitution into Eq. (73), with the result that 



ai(^)Pi(^) + a2(r)p2(r) 



ai[T)T- 



dpiir) 



a2[T)T 



•5^2 (r) _ p 



(75) 
(76) 



dr dr 

Solution of this pair of algebraic equations solves the problem of determining Gp. The Green 



function for u is then found by substitution into Eq. (|64|) to yield 

G^{k; r, r ) 



ir dGp{k; r, r') 



(77) 



k dr 

Without an explicit equation of state and viscosities it is not possible to be more specific. 



C. Linear equation of state 

To proceed further we now choose the equation of state P = Vs^e with a constant speed 
of sound Vs- This is a reasonable approximation to QCD at high temperature. The response 
functions Gp and can then be found in terms of elementary functions. 
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In the Langevin approach the drift matrix is a constant, and the evolution operator 
U can be found exphcitly in terms of the eigenvalues X± and eigenvectors il)± of D. With 



(78) 



we find (cf. Ref.[52]; 



X± =a ± f3 



13 =Va2 _y^2^2 



and 



Then the evolution matrix U can be determined as 



C/(A;;r,rO 



A+ -ik \ _ (rVr)^+ / A_ -ik 
-ikvJ^ -A_ / A+ - A_ I -ikvJ^ -X 



(79) 



(80) 



(81) 



The response functions are then given by Eqs. (72) as 



G,{k;T,r') 



Guj{k; T, t') = ik 



T 



T 



cosh(/31n(r/r')) + ( ) sinh (/31n(r/r')) 



cosh (/31n(r/r')) 



/3 



/3 



sinh(/31n(r/r')) 



(82) 



(83) 



Note that (3 is real if \k\ < (1 — fs^)/2t>s and is pure imaginary if \k\ > (1 — Vs'^)/2vs leading 
to exponential or oscillatory behavior, respectively. 



In the second method, the two solutions to Eq. ( 73 ) are found to be 



~p2{r) 



To \ "-/^ 



(84) 



The corresponding coefficient functions to construct G are 



a2[T 



2/3 V^o. 
P + a + k^ /r^^ °" 
2/3 Uo, 



(85) 



The response functions are identical to the ones given in Eqs. (82) and (83) 
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D. Singularities and the sound horizon 

It is instructive to analyze the singularities of Gp{$,;T,T') and G'(^(^; r, r'). As we shall 
see, they reflect propagation of sound along the z axis on top of the expanding medium. 
The propagation of sound waves in the transverse directions has been discussed in Ref. [8]. 

First, observe that Gp and Gi^ are meromorphic functions of k. The sole singularity is 
an essential singularity at infinity. As A; — t- oo, /3 — )■ ivsk whereas a remains a constant. 
Therefore, when |^| > fslog(r/r'), one can close the contour in the Fourier integral over k 
around either the upper or lower large semi-circle (depending on the sign of ^) to show that 

Gx{^;t,t') = when \^\>vMr/r'). (86) 

This means that there is a sound horizon which expands logarithmically with proper time r. 
This result is confirmed by the observation that, in the local rest frame the velocity of the 



front, rd^/dr, equals fg. Arguing similarly, or using Eq. (63) directly, one can show that 
the correlations do not spread beyond the sound horizon: 

Gxy{^; t, t') = when \^\ > 2v, ln(r/r') . (87) 

The singularities of Gxy{^\TiT') at the sound horizon can be analyzed by considering 
the large k asymptotics of its Fourier transform. For example, since for large k 

Gp{k;r,r') ^ - H-j sin [t;^^ ln(r/r')] , (88) 



we can use Eq. (62) to find that 

/°° dk ~ 
e^^^Gp{k-r,,T)Gp{-k-n, 
-oo 



4^ j ['^"(^ - + ^"(^ + Hn/r)) - 2^(0] (89) 

where ^ = S,2 — ^i- We see that the correlator is singular whenever there is a noise source 
event at earlier time r such that a sound cone originating from it goes through both points 
^1 and .^2 at time Tf. That source point is located midway between the points and .^2 for 



the first two delta-functions in Eq. (89). For the last delta-function, there are two source 
events located a distance ln(rf/r) away on either side of the coinciding points ^1 = ^2- 

The second derivative of the delta-function can be traced back to the derivative of the 
force term / with respect to ^. This derivative is a consequence of momentum conserva- 



tion. The second derivatives of the delta- function in Eq. (89) only represent the leading 
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singularities. There exist sub-leading singularities, such as step-functions, delta-functions 
and first derivatives of delta-functions. The prefactor (r/rf)^" describes the dilution of the 
fluctuation due to the expansion. 



E. The wake 



In a stationary medium, and neglecting viscous effects, sound propagation would be the 
only source of correlations. This is because the sound propagation would be non-dispersive 
with a linear relation between frequency and wavenumber. In the expanding medium we 
are considering, sound propagation also leads to a sound horizon but, unlike the stationary 



case, the dispersion relation is non-linear according to eq. (79). This leads to a wake behind 
the sound front, even without dissipation. 



Let us illustrate this using Eq. (82) to calculate the correlator Gpp{k;T{,T) 



\Gp{k; Tf, r)p. The Fourier transform Gpp{^] Tf, r) of this function is singular as discussed in 
the previous section. It is instructive to separate the most singular part of this function by 
expanding Gpp{k; Tf, r) in powers of and keeping all terms regular at k = 0. This leads 
to 

&;^^{k; Tf, r) = (aiP + b^) + [a^k^ + 62) cos {2v, ln(rf/r)A:) 

+ h^'-^^A sin (2^;, Hr,/T)k) , (90) 



k 

where the expansion coefficients a„ and 6„ are functions of ln(rf/r). The Fourier transform 



of Eq. (90), Gpp^{^;T{,T), is a sum of step functions and its derivatives with singularities 



located at ,^ = and C, = ±2fs ln(rf/r). The most singular term we have already written 



out in Eq. (89). The regular part 

G^pf ^ Gpp - (91) 

is a continuous function of ^, which is shown in Fig. [T| We see that the correlations spread 
along the rapidity axis from ^ = with time. In Appendix|A]we show that for long times this 
process resembles diffusion, and for asymptotically large Tf/r the function Gpp is given by a 
Gaussian. However, one should bear in mind that this diffusion is occurring in the absence of 
any dissipative effects since we have neglected viscous terms in the hydrodynamic equations 
at this stage. 
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0.30 




FIG. 1: The regular (continuous) part of the correlator G pp{^;T^^T) with v^^ = 1/3. The correlator 
is shown for ln(rf/r) = 2 (dashed line), ln(rf/r) = 4 (solid line), and ln(Tf/r) = 6 (dotted line). 
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FIG. 2: The singular part of the correlator Gpp{^;T{,T) with Vs^ = 1/3 and ln(rf/r) = 4 . The 
function is smeared by a Gaussian of variance cj^ = 0.1 in order to show the nature of the singu- 
larities. 



In order to display the singular part Gp^^ we convolute it with the Gaussian 



1 



-(?-r)2/2^2 



27r(T 



of a small width cr^ = 0.1. This simply replaces delta-functions with Gaussians and step 
functions with error functions. The result is shown in Fig. |2} 
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It is useful to note that the correlator Gpp{^; Tf, r) obeys the sum rule 

/oo 
d^Gppi^;n,r) = l, (92) 
-oo 

which is related to entropy conservation. Indeed, the time-independence of the integral in 



Eq. (92) is due to the zero mode in Eq. (78): X-{k = 0) = 0. The origin of the zero mode 



can be understood using the equation of motion (49). Expressing the fluctuation at A; = 
as 6p{k = 0, r) = const J 6s{C,, r) rdC,, we see that it is proportional to the fluctuation of the 
total entropy whose relaxation rate must vanish in the inviscid case. It is also interesting to 
note that at asymptotically large times Tf/r ^ 1 the sum rule is saturated by the regular 



part of Gpp since the singular part is suppressed by a factor (r/Ti) , see Eq. (89) and 
Appendix |X} 



F. Viscosity and taming of the singularities 



We now include the effects of viscosity on the space-time evolution of the system. In 
general, viscosity acts to smooth out gradients in temperature and flow velocity. Even if 
the viscosity is very small and its effects on the solutions to the equations of motion can 
mainly be neglected, the effect of viscosity on the correlation functions will still be important 
near the sound horizon singularities discussed in the previous two sections. The goal of this 
section is to demonstrate how the effect of viscosity smoothes out these singularities. The 



equations to be solved now are (53) and (54). For simplicity, and for illustrative purposes, 



we consider a constant value of u as well as a linear equation of state, vj^ = constant. 



As in the inviscid case, we combine Eqs. (53) and (54) into a two-component matrix 



Langevin equation. It is convenient to use a rescaled variable u {1 — v/Tt) instead of u. 



The drift operator now has the form D = Dq + Di with Dq being given by Eq. (68) and 
the viscous contribution (neglecting terms higher order in v/Tt) 



V 



-ik 



(93) 



To see the effect of viscosity more clearly it is convenient to rewrite the evolution operator 



in Eq. ( 71 ) as 



tj{k- r, r') = (ro/r)^« T exp [- £ ^ ^/(^")} (r' 



/To) 



Do 



(94) 
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where we defined the matrix Dj as 



Di{t") = (r7ro)^"r>i(r") (ro/r")^°. 



(95) 



Note that Di = 0{v). To leading order in v (in the exponent) we can use the Baker- 



Campbell-Hausdorff formula to remove the operation of time ordering in Eq. (94) which 



greatly simplifies calculations. The matrix TJ can be then calculated through a lengthy 
matrix algebra. 

In order to evaluate the effect of dissipation on the delta-function singularities in ^-space 
let us consider the limit /c — )■ oo. In this regime the matrix algebra simplifies and can be 



performed more easily. Keeping only the leading terms in Eq. (93) we obtain 



ik 
ikvJ 







(96) 



where cr^ denotes the third Pauli matrix. We next note that cr^ anticommutes with Dq in 



Eq. (96), which implies that Eq. (|95|) can be written as 



''7L[l--3(Vr")^"1 



(97) 



2T{t")t" 

Since the eigenvalues of Dq are pure imaginary (iifcfg), the second term in brackets in Eq. 



(97) is an oscillating function of Inr". Upon integration over t" in Eq. (94) the oscillating 



terms will be suppressed by a power of and the leading contribution, of order /c^, will 
be proportional to the unit matrix, an important simplification. 
The integral needed then is 

1 1 



^^'^ ' J^, t" T{t")t" 2a [T{t')t' T{t)t _ 



(98) 



where we used T = To(ro/r)"= and a defined in Eq. (79). Thus we can write 



U{k- r, t') ^ (ro/r)^« ^-uH{ry)kV2 /to)'''' = ^-^^-^-')^^ 1'^ Uo{k; r, r') (99) 

where LfQ{k]T,T') = (r'/r)-^° is the evolution matrix in the inviscid case. We see that the 
effect of viscosity is to dampen the oscillatory (or constant) behavior at large k as long as 
r 7^ t'. The overall result is that both Gp{k]T,T') and G^{k]T,T') are to be multiplied 



^ More precisely, the limit we consider is A; ^ 1 but vk^ ^ Tt. 
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by the Gaussian function Q-^H{T,T')k'^ /2 ^ TYie effect of viscosity is thus simply a diffusion in 
space-time rapidity whose proper-time dependence is controlled by the function H{t, t'). 
Next we turn to the second method for solving the same problem. In the limit of small 



viscosity, v <^ tT, Eqs. (53 54) in fc-space become 

dp 



T— + iku + f 
or 



T^ + {1- v,^)u + ik ( v,'p + / ) + 



-uo 



0. 



(100) 



(101) 



In the limit z/ = these reduce to exactly the same equations as studied earlier. The only 
new term is the last one in the second equation above on account of the fact that even 
though z//rT is assumed to be small, large enough values of k will make it important. Upon 
eliminating a) one arrives at a single second order differential equation. 



.2^ 



2 , ^^^\ 9p 2,2- 

2 - fs H — I r— +Vs k p 



tT J dr 



9f , /,2 



0, 



(102) 



Compared to Eq. ( 73 ) there is only one new term. 



To find the solutions to the homogeneous equation (/ = 0) it is convenient to change 
variables from r to a; = (tq/t)^". This leads to the second order differential equation 

d^p 



uk^ dp Vs^k"^ 

+ . o o p 



0. 



dx'^ 2aroTo dx Aa^x'^ 
The two independent solutions to this equation are 



(103) 



Pi 

P2 



xexp 
xexp 



4aroTo 
vk'^x 



/3/2a 



For uk'^/rT ^ 1, the arguments of both the Bessel functions and the exponential are small. 
Keeping the lowest order terms in both functions we obtain the same result in the inviscid 



vk'^x 



(104) 
(105) 



case as given by Eq. (84). When k is sufficiently large, vk'^/rT may not be small. In 
this regime, however, the index of the Bessel functions becomes large, (5/2a ~ 0{k), and 



the Bessel function can still be approximated^ as K^{z) 



z^. Therefore the 



^ This requires only that ^ /i, i.e. the argument of the Bessel function does not have to be small if the 
index is large. 
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solutions are approximately (the normalization doesn't matter) 



P2ir) 



a+l3 



exp 



(. 



exp 



V4arT(r) 
V4arT(r) 



(106) 



These are the same as the inviscid case Eqs. (84) with an additional exponential factor. 
The Green function Gp is constructed in exactly the same way as in the inviscid case, 



Eqs. (74 76) because the homogeneous term in / is unchanged. This now leads to 

a+/3 



ai[T 



a2[T 



1 
1 

2P 



P + a + k'^ + 







{- 


2t'T{t')_ 






uk^ ' 




fr' 


2T'T{r')_ 







exp 



uk' 



exp 



AaT'T{T'] 

uk^ 
' AaT'T{T') 



(107) 



The result for Gp is 



Gp{k]T, t') 



cosh (/3 ln(r/r')) 



/5 



sinh (/3 ln(r/r')) 



X exp 



vk^ 



4a 



t'T{t') rT{T))\ 

When z/ = it reduces to the inviscid case represented by Eq. (82). In the regime vk"^ /tT <^ 
1 the term proportional to v in front of the sinh is negligible and this result coincides with 



Eqs. (98| 99). The Eq. (108) is, however, somewhat more accurate, since it does not assume 
uk'^/rT <^ 1, which is reflected in the coefficient of the sinh. 



G. Example of other sources of smoothing of singularities 

The sound horizon singularities would be smeared if the delta-function correlator for 



the noise in Eq. (48) were replaced by a narrowly peaked regular function. Indeed, the 
origin of the noise is the fluctuation on a microscopic scale whose correlation length, small 
on the hydrodynamic scale, is non-zero. Here we shall consider the effect of the finite 
correlation length of the noise. Although the introduction of a finite correlation length is 
physically intuitive, it should be borne in mind that from the point of view of hydrodynamics 
it corresponds to inclusion of some, but not all, higher-order corrections in the systematic 
gradient expansion. We shall use this effect in the next section to estimate possible sensitivity 
of our results to higher-order hydrodynamic corrections. 
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The most microscopically sensible replacement for the delta-functions is one that is ex- 
ponential in time and in space. When distances are small the time and space intervals can 
be expressed in terms of Bjorken coordinates as At ^ Ar and Az ~ r A^. For simplicity, we 



keep the delta-function 6{t2 — ti) in Eq. (48) but replace 5(^2 — with an exponential 

^ ^ ^"'"'^'^^^ (109) 
r 2A(r) ' ^ ' 

where A(r) denotes the correlation length, which is a function of proper time. The net result 



is to multiply the right-hand side of Eq. (59) by 

^2 



rf + A2(ri)A;2 ' 

It is natural to assume that A(r) = cx/T{t). We can make a simple estimate of the constant 
c\ by assuming that A is given by the average interparticle distance at temperature T. For 
gluons plus three flavors of massless quarks the particle density is 



which gives c\ = n ^/^T = 0.637. The implications of this choice will be investigated in the 
next section. 



IV. PHENOMENOLOGY 



In order to make contact with experiment we need to consider how the fluctuations are 
frozen out. In general, hydrodynamic freeze-out occurs when the particles can no longer 
maintain local thermal equilibrium and therefore begin free-streaming. A schematic ap- 
proach to the freeze-out problem is represented by the Cooper-Frye formula [53] which 
describes the distribution of emitted particles as an integral over a freeze-out hypersurface 
Sf, usually chosen to coincide with a surface of constant temperature Tf (isothermal freeze- 
out) of the expanding fluid. In the Bjorken expansion scenario this surface is also a Tf = 
constant surface, but this equivalence holds only for averaged quantities. The fluctuations of 
temperature mean that the conditions Tf = constant and Tf = constant differ. In this paper 
we shall choose the simplest of these two conditions, isochronous freeze-out with Tf = con- 
stant, which has been used for the study of fluctuations in [51]. The alternative approach 
was pursued by Staig and Shuryak in their treatment of initial state fluctuations in the 
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transverse space [TO]. Both approaches offer only a schematic representation of freeze-out, 
but sufficient for our illustrative purposes. We leave the proper implementation of freeze-out 
to future studies. 



/ rfV^p^0(a-p)4/s(3^,p), (111) 



A. Freeze-out and rapidity smearing 

We start with the Cooper-Frye formula [S3] for the phase space distribution of emitted 
particles of a given species 

p dN, 
where 

/s(x,p) = (e(^-^=)/^±l)"' (112) 

denotes the local thermal distribution for the particle species with degeneracy ds-, and the 
step function ensures emission in the forward or outward direction. For the isochronous 
freeze-out in proper time the freeze-out hypersurface is given by r = Tf = constant. This 
means that the four-vector d^a^ has only a r-component in the Bjorken coordinates. Thus 
d^a ■ p = p'^TdC,d^x±, where the r-component of p^ is p'^ = m±cosh.{r] — C,) with m± = 
\/p\^ ml. Furthermore, 0{a^p'^) = 9{jf) = 1. For simplicity, we set /is = and neglect the 



quantum correction ±1 in the phase space distribution (112). For the purpose of obtaining 
numerical values we consider charged pions, dg = 2. We are interested in the number of 
particles of a given species per unit kinematic rapidity r] = tanh^^ (p^/pq). Integrating (111) 



over transverse area A and using {u'^,u^) = {cosh.u),T -^sinhw), {p'^,p^) = m_|_(cosh(?7 
^), sinh(r7 — 0)? obtain 



dN _^ r d'^p±A 



ds I fr,_\3 'Tfmj, / (i,^cosh(r7 — ,^) exp [— m_L cosh(?7 — ,^ — a;)/T] . (113) 



Fluctuations in dN/dt] are caused by fluctuations in temperature T (around Tf) and flow 
rapidity u. Expanding to linear order in fluctuations and integrating over d'^p± we find 

JdN\ d.AT,T^ f pv.^ + u tarih(r) -Q ^f, m„ \ 
H W j = 1^ J cosh^(.,-0 ^ V' T\ ''^') 

where T{x,y) is incomplete Gamma function and we used 

5T/T = v,^5s/s = v,^p. (115) 
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The rapidity correlator can then be written as 



drji dr]2 



dsATfT^\ 



(27r)2 



(116) 



X=p,LJ 

Y=p,ui 



where the smearing functions are given by 



Fpix) 



Fjx) 









cosh^(x) 






tanh(x) 


-r 




cosh^(x'; 







(117) 
(118) 



while Cxy{C,i — C,2] Tf) = {X{C,i, T{)Y{C,2, Tf) ) is the equal-time rapidity correlator defined in 



Eq. (61). It is convenient to use Fourier transforms of those functions in terms of which 



S—S— 
dT]i d-r]2 



where A?7 = ?7i — 772- 



I ^^..A, ^ Fx{-k)FY{k)Cxy{k-n), (119) 



X=p,u] 
Y=p,uj 



B. Normalization 



Since Cxy ~ 1/A (see Eq. (60)), it is convenient to divide by the event average of dN/drj 
given by Eq. (113) with a; = in order to remove the dependence on the transverse area A 
of the system. Integrating over d'^p± we get 

dN \ dsAr^T^ f dx f ^ ^0 



drj 



An' 



cosh^(x)n''^'"''^'' 



(120) 



Collecting all factors in front of the integrals in Eqs. (61), (119) and (120), one can write 
the normalized correlator as 



drji dr]2 



dN 
df] 



45d 



47r4iVeff(To) TfTf \T, 



7^2 



K{Ari) 



(121) 



where we defined the effective number of bosonic species such that s{T) = 2ti'^Ncq{T) T'^/45 
and To = T{tq). We also defined the dimensionless function K{Ari) in such a way that most 
of the dependence on To for the long-range tail is in the prefactor (using the observations at 



the end of Sec. [IIIEj). 

Assuming that reasonable values for the parameters are Vg^ = 1/3, T{ = 10 fm, Tf = 150 
MeV, To = 600 MeV, N^s = 47.5 (counting gluons and quarks) and u = l/3n, we plot the 
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normalized correlator (121) in Figs, [s] and |4j To determine K{^rj) we apply the freeze 



out (thermal) smearing described by Eq. (119) in both plots under the assumption that 



the observed particles are pions. To evaluate the effect of viscosity, we compare Fig. |3| 
which neglects viscosity, with Fig. |4] which includes viscous broadening as described by 



Eqs. (98) and (99). The factor in front of K(Ari) in Eq. (121) is approximately 1.1 x 10 



Finally, to estimate the effect of higher-order 



for our choice of parameters. Combining this with Fig. |4j, we conclude that Eq. (121) 
predicts correlations of the order of 10~^ 

3.0 
2.5 



<l 




FIG. 3: The correlation function K(Ar]) in the normalized correlator of dN/drj fluctuations in Eq. 



(126). Viscosity is not included. 



hydrodynamic corrections we consider the noise correlator with non-zero correlation length 
as discussed in Sec. IIIG Adding this effect on top of viscous broadening we obtain Fig. |5| 
The effect is visible but does not change the main features. 

The important conclusion is that the absolute magnitude of the correlation (outside of 
the 1^1 < 1 peak), for given freeze-out parameters, is proportional to the relative viscosity u 
and to a power of the initial temperature Tq. 



C. Contribution of initial state fluctuations 



Since the fluctuation equations (67) are linear and the noise is uncorrelated with initial 
conditions, the contribution of the initial state fluctuations to the correlator Eq. (61) is 
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FIG. 4: The correlation function K{Ar]) in the normahzed correlator of dN/di] fluctuations in Eq. 



(126). Viscosity is included. 

2.5 
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FIG. 5: The correlation function K(Ar]) in the normalized correlator of dN/di] fluctuations in Eq. 



(126). Viscosity is included, as well as the source correlator broadening discussed in Sec. IIIG 



additive and is given by 

CxY{k;Tty^' = J2 Uxx'{k;T,To)Cx'Y'{k;To)UYY'{-k;T,To) 



(122) 



X'Y' 



These correlations, unlike the purely hydrodynamic correlations discussed so far, depend also 
on the physics determining the initial-time correlator Cxy{C,', tq)- Calculation of Cxy{C,', tq) 
is beyond hydrodynamics; it could be done from the traditional Glauber approach or from 
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the color glass condensate (CGC) description of the initial state. However, once the initial 
correlator Cxvi.i^'TQ) is given, the subsequent evolution is governed by hydrodynamics ac- 



cording to Eq. (122). Since hydrodynamic evolution is the main subject of this paper, we 
shall assume a generic form of the initial correlator Cxy{^'-, ^o), leaving its calculation beyond 
the scope of the paper. It is reasonable to assume that this correlator is local, meaning that 
Cxvik', To) is a polynomial in k. We shall also assume, for simplicity, the following matrix 
form for it which obeys the basic symmetry properties of the correlator 

1 -tk 



CxY{k;To) = {X{k;To)Y{-k;To)) 



Co 

A 



ik fc^ 



(123) 



XY 



where the factor 1/A is due to the locality of the correlator in the transverse space (Eq. 



(60)) and where we defined the dimensionful coefficient cq which parameterizes the absolute 



strength of the correlator. Substitution into Eq. (122) results in 

Co. 



CxY{k;Ty 



A 



GxY{k;T,To) 



(124) 



where we used the definition of Gxy in Eqs. (62) and (72). Up to the constant factor cq, 
the Fourier transform of this correlator has been already discussed and its matrix element 



G pp has been plotted in Sec. IIII E 



The contribution of such initial state fiuctuations to the two-particle rapidity correlation. 



similarly to Eq. (|119), is given by 

i 



S—S— 
drji dr]2 



47r2 ) 



dk 
— f 
27r 



X=p.u) 

Y=p',LJ 



Fx{-k)FY{k)GxY{k; Tf, ro) . (125) 



Collecting the factors in front of the integrals in Eqs. (125) and (120) we can write the 



contribution of initial state fiuctuations to the normalized correlator as 

dN\~^ _ d^TiTf 
drj / ° A-k"^ 



^(M ' '''' 
dr]i dr]2 



Co 



ir™(Ar/) . 



(126) 



Here we defined the function which we plot in Fig. p\ for the same choice of the 



parameters as in Sec. IV B The coefficient in front of K"^\At]) in Eq. (125) is 8.5 x 
10~^(co/l GeV^) for that choice of parameters. Assuming that CGC initial conditions give 
rise to cq of the order characteristic saturation scale cq ~ QLt ~ ^ GeV^ we conclude that 



such initial state fiuctuations produce correlations of similar magnitude to those due to 
purely hydrodynamic fiuctuations. Of course, a more detailed analysis of initial conditions 
is needed before a quantitative comparison with experiment can be made. 
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FIG. 6: The correlation function K^^^{Arj) in the normahzed correlator in Eq. (126). Viscosity is 
included. 



V. SUMMARY AND CONCLUSIONS 



In this paper we explored the contribution of local hydrodynamic fluctuations to the 
event-by-event fluctuations of particles emitted from relativistic heavy ion collisions. Un- 
like the contribution of initial state fluctuations, which have been the main focus of the 
studies so far and whose magnitude is determined by quantum pre-equilibrium dynamics, 
the magnitude of the fluctuations we discussed here is directly related to the hydrodynamic 
properties of the locally equilibrated matter. In the framework of relativistic viscous hydro- 
dynamics, owing to the fluctuation-dissipation theorem, the amphtude of these fluctuations 
is governed by the viscosities. This offers a possibility to measure, or constrain, the viscosity 
of the strongly coupled quark-gluon plasma independent from the traditional analysis of 
elliptic flow. 

We observed two remarkable features of the fluctuation correlator. The first is that 
the correlations spread in rapidity space logarithmically with Bjorken proper time, with 
velocity determined by the speed of the sound in a static medium. This behavior is similar 
to the circles observed in [9] for correlations in transverse space induced by initial state 
fluctuations. The second is that we find the correlations are not limited to the sound cone 
but are accompanied by a wake behind the sound front which can be traced to the non- 
linearity of the sound mode dispersion relation in the medium. This diffusion- like wake is a 
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non-dissipative process. 

At the lowest order in the gradient expansion, the source of the fluctuations is white 
noise, which naively leads to singularities in the correlation functions of the density and 
of the flow velocity. The singular behavior is tamed, in part, by the viscous terms in the 
hydrodynamic equations, as well as by the thermal smearing when one calculates final-state 
particle distributions. For completeness, we also considered the effect of replacing the white 
noise by colored noise with a thermal correlation function. 

We explored the phenomenological consequences of hydrodynamic fluctuations in the 
idealized scenario of boost-invariant longitudinal flow with a homogeneous transverse profile. 
In the Cooper- Frye approach to freeze-out, the correlation function of the particle yield 
dN / drj as a function of the kinematic rapidity difference A?7 is obtained from the temperature 
and flow velocity correlation functions in Bjorken space rapidity by a thermal smearing. 

Two features of the particle number correlation function deserve special mention. One 
is a strong peak at Ar] = 0, which receives contributions from hydrodynamic fluctuations 
during the entire course of the expansion. Its height and width is influenced by both the 
thermal smearing at freeze-out and the viscous smearing during the hydrodynamic phase 
(see also [55J). 

The other noteworthy feature is a broad structure at larger rapidity differences, extend- 
ing up to the sound horizon Ar/max = 2i;s ln(rf/ro), which is caused by the hydrodynamic 
propagation of the noise followed by a slower diffusive wake generated at early times until 
thermal freeze-out. Depending on the precise values of the sound velocity and the start and 
end of the hydrodynamic phase, this implies that particle number correlations extending 
over significantly more than one unit of rapidity can be generated during the hydrodynamic 
phase. Correlations over large rapidity intervals have been observed in heavy-ion experi- 
ments [301 EI] cind have been subject to numerous theoretical studies [321439] . It would be 
interesting to investigate to what extent the hydrodynamic fluctuations contribute to this 
phenomenon. 

Let us emphasize again the main point of this paper. It has been clear for some time that 
the profile of the particle number correlations depends on the values of the shear (and bulk) 
viscosity and of the sound velocity. As others have already argued [H [28] this is true for 
azimuthal correlations generated by fluctuations in the hydrodynamic initial conditions over 
the transverse plane. Our results conflrm this for the longitudinal space correlations. What 
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distinguishes hydrodynamic correlations induced by local thermal fluctuations is that their 
absolute magnitude is also determined by hydrodynamic properties of the medium. Thus 
our findings amplify the opportunities offered by measured particle number fluctuations to 
constrain the fluid dynamical properties of the hot matter created in relativistic heavy ion 
collisions. 
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Appendix A: Long-time limit 



Here we shall determine the behavior of the correlator Gpp{^; Xf, r) in the limit of asymp- 
totically large time separation rf/r ^ 1. We shall be able to obtain a closed analytical 
expression for the correlator in this limit. Unfortunately, for realistic values of Tf/r relevant 
for heavy-ion collisions this expression is still a poor approximation. However it is still use- 
ful, as analytic solutions often are, in demonstrating conceptually important features of the 
fluctuation correlator. 

In the limit we consider the behavior of the correlator is determined by the modes with 



the slowest rate given by eigenvalues (79). We observe that the slowest mode is X-{k) at 
small k. This mode relaxes arbitrarily slowly as A; — )■ like 



A. 



(Al) 



where 7^ = 1/ a/1 — v^, while the mode Xj^ = \ — v^,^ + 0{kP'). Thus in the limit r/r' ^ 1 
the longest lingering modes are smooth modes corresponding to the eigenvalue A_ given by 



in Eq. (80). For such a mode 



UJ = ^ -ikvl^ilp . 



(A2) 



The matrix U in this limit becomes a projector on this mode, as can be seen directly from 



Eq. (81 ). Instead of proceeding from there, we can also simply substitute Eq. (A2) into Eq. 



(100) and obtain a single equation for p 



r^ + X_p + f = 

OT 



(A3) 



which is easily solved in the form of Eq. (57) with the Green's function Gp{k; r, r') given by 



Gp{k;T,T') = {T'/Ty 



(A4) 



Using equation ( 62 ) we find 



Gpp{k;Tf,T) = (r/Tf) 



(A5) 



where we also used (|Al|). Fourier transforming from A; to ^ we find 

1/2 1 



Gppi^; Tf, r) = (87rfs^7,^ ln(rf/r)) exp 



e 



8fs2^2 ln(rf/7 



(A6) 
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We see that the fluctuation correlator describes diffusion in the Bjorken coordinate ^ with 
ln(rf/r) playing the role of time. It is interesting that this diffusion-like process is not 
associated with dissipation. 

Another question which can be asked is what happens to the singularities we discussed 



in Sec. HID in the long-time limit. This can be easily seen by comparing Eq. (A6) with Eq. 



(89). One can see that the strength of the singularities decreases exponentially with ln(rf/r) 



as (r/rf 



, while the contribution of the smooth modes is roughly time- independent. In 



particular, the sum rule (92) is completely saturated by the Gaussian (A6) at late times. 



Numerically, at finite ln(rf/r) = 4 the regular part G^^^ in Fig. [Tj contributes 1.12, with the 
singular part Gp™^ accounting for —0.12. 



Appendix B: Azimuthal correlations and power spectrum 

Our example application of the theory of hydrodynamic fluctuations, the one-dimensional 
boost invariant Bjorken flow, allowed us to study correlations of the hydrodynamic quantities 
in the longitudinal direction. There are several natural extensions. Here we shall make a 
few comments concerning the study of hydrodynamic correlations in the transverse plane. 
The natural quantity to consider in this case are azimuthal correlations among the number 
of emitted particles event by event. 

In contrast to initial state fluctuations, which are intimately related to geometric aspects 
of the nuclear collision, hydrodynamical fluctuations are caused by local noise and are thus 
unaffected by the global geometry. This suggests that it makes little sense to look for corre- 
lations between hydrodynamical fluctuations and global event properties, such as the event 
plane deflned by the impact parameter or the elliptic flow pattern. Instead, a more promis- 
ing analysis will follow the procedure used to determine the power spectrum of fluctuations 
in the cosmic microwave background radiation (CMBR) [56]. The general argument for 
the application of this approach to heavy ion collisions has been proposed by Mishra et al. 
[57] [58] for the elliptic flow velocity V2- 

The experimental determination of the power spectrum of azimuthal fluctuations relies on 
the measurement of two-particle correlations in the flnal state. This observable was originally 
suggested as a method of measuring collective flow anisotropics that does not require the 
determination of the reaction plane |59H6T] . However, as Mishra, et al. emphasized, the 
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power spectrum |f„p of the azimuthal anisotropy of the distribution of emitted particles 
not only picks up the collective flow anisotropy, but also event-by-event fluctuations of 
the emission pattern, including density fluctuations and fluctuations of the collective flow 
velocity. 

Within a chosen rapidity window, which will depend on the masses of the particles in- 
volved, one can represent the distribution in some observable O as 

O(0) = O + 5^o„e^"^ (Bl) 

where denotes the azimuthal angle. The window is fixed in the laboratory frame of 
reference, and its orientation does not vary from one collision to the next. For the CMBR 
the averaging is performed over points in the sky. In the case of relativistic heavy ion 
collisions the averaging is done over a large set of events in which the orientation reaction 
plane and the hydrodynamic noise change randomly from one event to another. Thus 

(o„) = 0, {or^ol,) = \0l5^y , (B2) 

where 0„ characterizes the magnitude (power) of the rfi^ angular Fourier component of the 
fluctuations. The dependence of 0„ on n may reveal hydrodynamic and thermodynamic 
properties of the expanding medium. For example, as has been already observed in the 
studies of correlations induced by initial state fluctuations [9j, viscosity suppresses higher 
harmonics n. It would be also interesting to consider the effect of the QCD critical point 
(see Ref. [62] for a review) on the power spectrum. In this case the natural choice of the 
variable O would be the baryon density, or its experimental proxy - net proton density. The 
magnitude of fluctuations increases near the critical point [HSl El] , and the fluctuations can 
become highly non-Gaussian [651 ES]. The increase of fluctuations could be observed in the 
power spectrum. In addition, the increase in the correlation length may cause the power 
spectrum to shift its weight towards smaller values of n. Further quantitative investigation 
is needed, of course, to determine whether these effects have an observable magnitude. We 
leave this to future work. 



